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We report results of a high-resolution numerical-relativity simulation for the merger of black hole- 
magnetized neutron star binaries on Japanese supercomputer “K”. We focus on a binary that is 
subject to tidal disruption and subsequent formation of a massive accretion torus. We find the 
launch of thermally driven torus wind, subsequent formation of a funnel wall above the torus and 
a magnetosphere with collimated poloidal magnetic field, and high Blandford-Znajek luminosity. 

We show for the first time this picture in a self-consistent simulation. The turbulence-like motion 
induced by the non-axisymmetric magnetorotational instability as well as the Kelvin-Helmholtz 
instability inside the accretion torus works as an agent to drive the mass accretion and converts the 
accretion energy to thermal energy, which results in the generation of a strong wind. By an in-depth 
resolution study, we reveal that high resolution is essential to draw such a picture. We also discuss 
the implication for the r-process nucleosynthesis, the radioactively-powered transient emission, and 
short gamma-ray bursts. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 


I. INTRODUCTION 

The merger of a black hole (BH) and a neutron star 
(NS) is one of the most promising sources of gravita¬ 
tional waves m- It could be also one of the strongest 
high-energy phenomena in the universe, if the NS is 
tidally disrupted by the companion BH before the on¬ 
set of the merger. Previous numerical-relativity simu¬ 
lations [IHH] have shown that after tidal disruption, a 
system composed of a BH and a torus is formed. This 
system could be the central engine of short gamma- 
ray bursts (sGRBs) [IS]. It is also found that dur¬ 
ing tidal disruption, an appreciable amount of mass is 
ejected laisKioi [TT] . Such ejecta, which should be 
neutron-rich, will subsequently produce heavy elements 
via r-process nucleosynthesis HZ]. The produced un¬ 
stable heavy elements will subsequently decay, heat-up 
the ejecta, and shine [Hj. Such an electromagnetic sig¬ 
nal could be an electromagnetic counterpart to detected 
gravitational waves [MlZj- 

These facts motivate us to perform physically reliable 
numerical-relativity simulation of BH-NS mergers. Here, 
the presence of strong magnetic fields is one of the most 
characteristic properties of NSs [5S]. However, the role of 
the magnetic fields in their merger process is still poorly 
known. 

BH-NS mergers can be subject to tidal disruption, for a 
broad range of the NS compactness, the mass ratio of BH 
to NS, and BH spin HHHj- An accretion torus, expected 
to be formed around the remnant BH, is subject to the 
magnetorotational instability (MRI) [29|, and thus, the 
magnetic field will be amplified. Previously, it was dif¬ 


ficult to perform high-resolution simulations to resolve 
the fastest growing mode of the MRI because of lim¬ 
ited computational resources, although several prelimi¬ 
nary numerical-relativity simulations have been carried 

out [Mai]. 

In this paper, we report the results of our latest gen¬ 
eral relativistic magnetohydrodynamics (GRMHD) sim¬ 
ulation for the BH-NS merger performed on Japanese 
supercomputer K. The highest-resolution simulation per¬ 
formed so far together with an in-depth resolution study 
was done. 


II. METHOD, INITIAL MODELS, AND GRID 
SETUR 

Einstein’s equation is solved in a puncture-Baumgarte- 
Shapiro-Shibata-Nakamura formalism together with 
fourth-order finite differencing [35U88] and the GRMHD 
equations are solved by a high-resolution shock captur¬ 
ing scheme. The simulations are performed using a 
fixed-mesh refinement (FMR) algorithm in which each 
refinement level labeled by i covers the cubic domain 
of G [—A^Ax(j),fVAx(j)] with Ax(i) being the grid 
spacing of level i. Ax(q = 2Ax(i_|_i) and i = 1,2 • • • , 
and fmax — I (see Refs, po] 00] for details). Typically, we 
set imax = 10 and the finest grid domain is a (123 km)^ 
cube. To examine how the result depends on the grid 
resolution, we change Aximax = 120m, 160m, 202m, and 
270m, respectively. We show the grid set-up in our sim¬ 
ulations in Table m In the highest-resolution run, we use 
32, 768 CPUs. 
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As initial data, we prepare a BH-NS binary in quasi¬ 
equilibrium using the method of Ref. [41]. We model 
the NS by the Akmal-Pandhalipande-Ravenhall EOS jUj 
which is compatible with a maximum neutron star 
mass > 2Mq as required by current observational con¬ 
straints |43l |44|. We set the NS mass, the mass ratio 
of BH to NS, and the dimensionless aligned spin of BH, 
X, to be I. 35 M 0 , 4, and 0.75, respectively. The initial 
orbital angular velocity is Gmofl/c^ = 0.056, where G 
is the gravitational constant, mo is the sum of BH and 
NS gravitational mass in isolation, and c is the speed of 
light. With these parameters, a massive accretion torus 
is formed after tidal disruption |5| j^Sj . 

The initial magnetic-field configuration is given in 
terms of the vector potential as |15| 

= {-(y - yNs)Sj +{x- Xns)^]') max(P - 0 )^, 

where a:NS and ?/ns denote the coordinate center of the 
NS, P is the pressure, Pc = P{p = O.Odpmax), and 
j = x,y, and z. Pmax is the maximum rest-mass density 
and we set such that the initial maximum magnetic- 
field strength is 10^® G. Even if we start a simulation 
with this ad hoc localized seed magnetic field, the re¬ 
sulting torus surrounding the remnant BH becomes in a 
turbulence-like state and a global magnetic field is nat¬ 
urally formed eventually. The magnetic-field strength 
is chosen so that the wavelength of the fastest grow¬ 
ing mode of the non-axisymmetric MRI is larger than 
lOAxiinax- We note that the resulting turbulence-like 
state should not depend critically on the initial magnetic- 
field configuration and strength as long as the MRI is re¬ 
solved HD (see also Ref. m for the discussion of the de¬ 
pendence of the saturation amplitude on the initial field 
configuration). The EOS is parametrized by a piecewise 
polytrope [35] and a gamma-law EOS is added during the 
simulation to capture shock heating effects with thermal 
gamma of 1.8. Other choice of thermal gamma could 
affect the amount of the torus wind. 


III. RESULTS. 


Eigurejl] shows snapshots of the rest-mass density pro¬ 
file with the magnetic-field lines. Before swallowed by 
the BH, the NS is tidally disrupted. A part of NS mat¬ 
ter subsequently forms an accretion torus around the 
remnant BH with mass of « O.ISMq at « 10 ms af¬ 
ter the tidal disruption. The Kelvin-Helmholtz instabil¬ 
ity develops in the contact interfaces of the wound spi¬ 
ral arm because of the presence of shear motion shown 
in Fig. [T}i and the vortices are generated subsequently 
enhancing the magnetic-field energy. In addition, the 
non-axisymmetric MRI activates amplification of the 
magnetic-field strength [49] (Fig. The mass accre¬ 
tion is enhanced by turbulence-like motion which is gen¬ 
erated by these MHD instabilities as well as by the grav¬ 
itational torque exerted by the non-axisymmetric struc¬ 


ture of the accretion torus (see Fig.[^ and visualization 
in Ref. [50]1. 

Figure plots the ejecta mass, torus mass, mass ac¬ 
cretion rate onto the BH, and ratio of the magnetic-field 
energy, E-q, to internal energy, E-mt, as functions of time. 
We define the merger time tmrg to be the time at which 
the gravitational-wave amplitude becomes maximal. The 
ejecta is defined to be fluid elements which reside outside 
the BH and have Ut < —1 where Ut is the lower time 
component of the four velocity. This criterion means 
that the fluid elements are gravitationally unbounded. 
The primary mass-ejection mechanism is tidal torque ex¬ 
erted during tidal disruption and this dynamical ejecta 
is seen for 0 ms < t — tmrg ^ 10 ms in Fig. During 
this early phase, NS matter of ~ O.OIMq is ejected ap¬ 
proximately along the orbital plane. This result agrees 
with that found in Refs. [Hi, which demonstrates that 
our findings will remain robust also for a larger initial 
separation. 

After this primary phase, a new ejecta component ap¬ 
pears. In the highest-resolution run, the accumulated 
accretion mass onto the BH for 10 ms < t — tmrg ^ 30 
ms is « O.O 3 M 0 , while the torus mass decreases by 
« O.O 6 M 0 over the same time. This implies that a signif¬ 
icant amount of the torus mass is ejected by a torus wind. 
The launch time and amount of material ejected by the 
wind depend strongly on the grid resolution: The higher- 
resolution runs result in an earlier launch time and larger 
amount of the ejecta. The mass accretion rate onto the 
BH also depends on the grid resolution: For higher res¬ 
olutions, it is smaller. The reason for these facts will be 
described later. 

The bottom panel of Fig. shows that irrespective 
of the grid resolution, the magnetic-field energy is ex¬ 
ponentially enhanced and eventually saturated: Eb is 
typically 0.1% of Ei^t- The growth rates of Eb for 
10 ms < t — fmrg ^ 20 ms correspond to 7-8% of the 
orbital angular velocity. This growth rate agrees approxi¬ 
mately with that of the non-axisymmetric MRI predicted 
by the linear perturbation analysis |49|. 

To see that our grid setting is sufficient to resolve 
the fastest growing mode of the non-axisymmetric MRI, 
Fig.j^plots a snapshot of wavelength of the fastest grow¬ 
ing mode of the non-axisymmetric MRI on a meridional 
(x-z) plane at < — tmrg « 15.0 ms. We estimate the wave¬ 
length by 


Amri,m = 


Av) 


2tt 


where 6 (,^), p, h, and H are an azimuthal component of 
the magnetic field measured in the fluid rest frame, the 
rest-mass density, the specific enthalpy, and the angular 
velocity, respectively. The wavelength is longer than « 3 
km in a large portion of the region and this indicates 
the fastest growing mode is covered by more than ten 
grid points even in the lowest-resolution run. The right 
panel of Fig. clearly shows it. Therefore, turbulence¬ 
like motion produced by the MRI, which is resolved in 
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our numerical simulation, plays an important role in the 
mass ejection. We discuss this point later. 

In the presence of neutrino radiation, the growth rate 
of the fastest growing MRI mode could be significantly 
suppressed once neutrino viscosity and drag turn on . 
According to Ref. [53], the neutrino viscosity z/, the neu¬ 
trino mean free path Ai/_mfp, and the wavelength of the 
fastest growing MRI mode AvIs-mri in the viscous regime 
are : 

V = 1.2x 10^^ TIq s“^, 

Ai/—mfp — 10 Pi 2 ^10 cm, 

Avis-MRI = 2.4 X 10^ cm, 

where Tip = T/lOMeV, pi 2 = p/I 0 ^^gcm“^, Ila = 
n/10^s“^, and z^i 2 = z^/10^^cm^s“^. Note that these 
estimates would depend on the structure of the accretion 
torus [54l[55] . 

Assuming the gas, photons, and relativistic electron 
and positrons contribute to the specific thermal energy, 
we evaluate the temperature from the thermal compo¬ 
nent of the specific internal energy |^. We utilize the 
data of the density, angular velocity, and the specific 
internal energy along a;-axis on the equatorial plane at 
t — tmrg ~ 10 ms. 

Figure ffl plots the radial profile of Ay_i„fp and 
Avis-MRI- In the entire region, A^-mip is always longer 
than Avis-MRI, which implies that the effect of neutri¬ 
nos on the MRI is not described by the viscosity, but 
by the neutrino drag. The neutrino drag is character¬ 
ized by the damping rate F of the velocity fluctuation 
due to the momentum transport. According to Ref. |S3], 
F = 6 X lO^TfgS”^. The radial profile of F is shown in 
Fig.i Because F ^ 11 in the entire region, the growth 
rate of the MRI is not affected by the neutrino drag. Al¬ 
though there is an ambiguity in terms of the accretion 
torus structure and the temperature estimation, we con¬ 
clude that the MRI growth rate is not significantly differ¬ 
ent from that of the ideal MHD in this BH-NS merger. 
The MRI will grow exponentially even if we assume a 
weak magnetic-field strength of ^ 10^^ G. 

Figure plots snapshots of the rest-mass density, 
plasma /3 (the ratio of matter pressure to magnetic pres¬ 
sure) , thermal component of specific internal energy, and 
sum of the Maxwell and Reynolds stress on the x-z plane 
at < — tmrg ~ 50.6 ms for the highest-resolution run. We 
also plot contours of Ut. Here, W = ln(— m*) is an effec¬ 
tive potential of a test particle moving in stationary and 
axisymmetric BH spacetime m- The shape of curves 
in the vicinity of the rotational axis with ut = —1 is 
approximately parabolic. In the Newtonian limit, W is 
reduced to —GM'byi/{R^+z^Y/'^+P/2R^ where Mbh, 
and I are the BH mass, cylindrical coordinate, and spe¬ 
cific angular momentum, respectively. We assume con¬ 
stant specific angular momentum for simplicity. In this 
case, for a given value of I the contour of m* = —1 be¬ 
comes parabolic. If the specific angular momentum for 
fluid elements is sufficiently enhanced or fluid elements 


are pushed to high latitude {z > R) by thermal pressure, 
they could have ut < — 1 {W > 0) [tS] (ST] [58] (see the 
discussion in the next paragraph). Because there is no 
matter in the region above the torus, the wind, once it 
is launched, expands in the widely spread radial direc¬ 
tion by contrast with the tidally induced ejecta. Subse¬ 
quently, a funnel wall is formed. 

The point to be clarified is how fluid elements are in¬ 
jected into the region with ut < —1 {W > 0). We find 
that j3 at the launch time of the wind is much greater 
than unity near the torus; pure magnetic pressure would 
not be the main agent of the injection. The bottom- 
left panel of Fig. indicates that there is a hot region 
in the vicinity of the BH, which produces a steep pres¬ 
sure gradient. Due to this gradient, the fluid elements 
are accelerated radially and become unbound once they 
reach the region with ut < —I {W > 0)(see Ref. [52] for 
essentially the same discussion). 

To explore the mechanism to enhance the thermal 
pressure, we analyze the specific kinetic-energy spectrum 
E{k), which is calculated by 1/2 / {Sv^(x + 

r)5v^{x))<PrdQ,k where fc is a wavenumber vector, k = 
\k\, dVlk is the volume element in a spherical shell between 
k and k -|- dk, and R is a cubic region of x[km] e [50, 70], 
?/[km] e [—10,10], z[km] S [—10,10]. We choose x as the 
center of the cube and r is the position vector from the 
center. The fluctuation of the velocity du* is u® — (u*) 
inside the chosen cube. (•) denotes the time average for 
the duration 10 ms < t — < 20 ms. 

Figure]^ plots this energy spectrum for all runs. This 
shows that (i) the spectrum is extended to smaller scales 
in the higher-resolution runs, indicating that the en¬ 
ergy is injected at a small scale at which the MRI de¬ 
velops and (ii) the spectrum amplitude is higher for 
the higher-resolution runs. Numerically resolvable wave¬ 
length of the fastest growing mode of the MRI becomes 
shorter with increasing the grid resolution. Therefore, 
the wavenumber for which the MRI develops is given 
by k/2TT = l/lOAximax ~ 0.83 km“^ for the highest- 
resolution run, « 0.62 km~^ for the middle-resolution 
run, « 0.50 km~^ for the normal-resolution run, and 
« 0.33 km“^ for the lowest-resolution run, respectively. 
This is expected to be approximately equal to the en¬ 
ergy injection scale in E(k). Assuming that a turbulent 
state is realized, the specific energy dissipation rate is 
~ Sv^/ledd where 4dd is the scale of the vortices [H2]. Be¬ 
cause the turbulent energy is proportional to 6v'^, Fig. 
suggests that Sv becomes higher in the higher-resolution 
runs for a given scale ledd- This indicates that the effec¬ 
tive turbulent viscosity, leddSv, increases with increasing 
resolution. Due to the realistic high viscosity achieved by 
this turbulence-like motion, the mass accretion inside the 
accretion torus is enhanced in the higher-resolution runs 
and the mass accretion energy is efficiently converted to 
thermal energy in the vicinity of the BH |63II65| . 

Figure]^ plots snapshots of the thermal component of 
specific internal energy on the x-z plane at t—tmrg ~ 50.6 
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ms for Ax = 160 m and 270 m. Comparing to the 
bottom-left panel of Fig. the thermal component of 
specific internal energy in the vicinity of the inner edge 
of torus and the torus interior is larger than that in the 
rest of torus in the higher resolution run. This trend 
indicates that efficient thermalization of the mass accre¬ 
tion energy due to the effective turbulent viscosity would 
be realized as discussed above. Consequently, the torus- 
wind power is enhanced while the mass accretion onto 
the BH is suppressed [68H65] . which can be seen clearly 
in Fig. The bottom-right panel of Fig. [^indeed shows 
that the energy is transferred outward. In Figs. and 
the amount of the ejecta mass and the spectrum of the 
matter flow do not exhibit the convergence. However, 
with the improvement of the grid resolution, the total 
amount of the torus wind mass increases and this indi¬ 
cates that our highest-resolution results would show the 
lower bound of the total wind mass. 


The high BH spin with a small horizon radius, which 
is necessary for tidal disruption of realistic BH-NS bina¬ 
ries [a [MU], prevents the fluid elements from being ac¬ 
creted on the remnant BH. Then, the fluid elements tend 
to stay in the vicinity of the BH and the pressure gradi¬ 
ent is enhanced EOj. This is also the key to push the fluid 
elements outward |52j . In our simulation, the BH is spun 
up to X « 0.9, which implies the radius of the inner most 
stable circular orbit is so small « 2.32GMbh/c^ (SH) that 
an efficient draining by the BH is prohibited. Hence, the 
accretion onto the BH is suppressed (BO] [ 6 T| (see also the 
panel (c) of Fig. [^. The amount of the torus wind mass 
is « O.O 6 M 0 in the highest-resolution run, which corre¬ 
sponds to about 50 % of the torus mass at t — fmrg ~ 10 
ms. 


As discussed above, the pressure gradient in the vicin¬ 
ity of the BH accelerates the outflow and this results in 
the formation of coherent poloidal magnetic helds be¬ 
cause the magnetic-field lines are frozen into fluid ele¬ 
ments. A low-/3 region is formed along the funnel wall 
in the wind phase. Subsequently, the magnetic pressure 
pushes matter and magnetic-field lines to the polar re¬ 
gion because there is only dilute matter in this region 
at the wind launch. This results in the formation of a 
BH magnetosphere. The top-right panel of Fig. [^indeed 
shows that a region with /3 ~ 10“^ is formed around z- 
axis. In the presence of a BH magnetosphere composed of 
a coherent poloidal magnetic field, the Blandford-Znajek 
(BZ) mechanism efficiently works for the outgoing 
Poynting flux generation. Figure [^ plots time evolution 
of the outgoing Poynting flux estimated on an apparent 
horizon. This figure shows that the Poynting flux is sig¬ 
nificantly enhanced after the wind launch because of the 
coherent poloidal magnetic-field formation.The Poynting 
flux is as high as « 2 x 10 "*® erg/s in the end of the 
highest-resolution run. 


TABLE I. Grid set up for four different grid-resolution runs. 
*max: The number of the refinement levels. : The 

grid spacing of the finest refinement level. N: The grid num¬ 
ber in one positive direction. 


Model 

^max 


[m] N 

high 

10 

120 

514 

middle 

10 

160 

378 

normal 

8 

202 

306 

low 

8 

270 

224 


IV. SUMMARY AND DISCUSSION. 

We performed high-resolution GRMHD simulations 
of a BH-NS merger on the supercomputer K. We self- 
consistently show a series of the processes composed of 
tidal disruption of the NS, the accretion torus forma¬ 
tion, the magnetic-field amplification due to the non- 
axisymmetric MRI, thermally driven torus wind, sub¬ 
sequent formation of the funnel wall and BH magneto¬ 
sphere, and the high BZ luminosity. 

A resolution study revealed that turbulence-like mo¬ 
tion works as the agent to drive the mass accretion and 
convert kinetic energy to thermal energy resulting in the 
generation of a strong wind. To show this phenomenon, 
sufficiently high-resolution simulations are essential. Af¬ 
ter the launch of the torus wind, a funnel wall and mag¬ 
netosphere with collimated poloidal magnetic fields are 
naturally formed. 

The torus wind and subsequent funnel plus magneto¬ 
sphere formation have the following implications. First, 
the formed magnetosphere could help launching a jet by 
the BZ mechanism. The high outgoing Poynting flux 
found in our simulation could be the main engine for 
sGRBs [BTj. Also, the jet could be collimated naturally 
by the pressure exerted by the funnel wall once the wind 
is launched. 

The torus wind could contribute significantly to r- 
process nucleosynthesis of heavy elements in BH-NS 
mergers. The dynamical ejecta would be neutron-rich 
and have a low value of electron fraction Ye- By con¬ 
trast, the torus wind component is expected to have 
a higher value of due to weak interactions because 
it has high temperature by shock-heating | 68 j . A mix¬ 
ture of the dynamical and wind components could be a 
key to reproduce the solar abundance pattern of the r- 
process heavy elements. Note that it was suggested that 
viscosity-driven and neutrino-driven winds from a torus 
around the BH could reproduce the solar abundances for 
mass-number greater than 90 (BS) (see also Refs. iMiiini 
for the NS-NS case). 

Finally, we comment on the kilonova/macronova 
(radioactively-powered electromagnetic emission) 
model |18j . The amount of the torus wind component in 
the highest-resolution run is as high as « 0 . 06 Mq which 
is much larger than that of the dynamical component 
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FIG. 1. Snapshots of the rest-mass density profile with the magnetic-field lines (a) just after tidal disruption , (b) at an early 
phase of accretion torus, and (c) in the final phase. The iso-surfaces for 10^^, 10^°, and 10® g/cm® are denoted by yellow, green, 
and blue color. The magnetic-field lines are shown by the white curves. In the middle panel, the iso-surfaces are drawn for the 
three-quarter region. 



FIG. 2. Time evolution of (a) ejecta mass, (b) torus mass, 
(c) mass accretion rate onto the BH, and (d) magnetic-field 
energy divided by internal energy. In panel (c), the data for 
t — tmrg JS 20 ms for Ax = 202 m run is not plotted. 
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FIG. 3. Profiles of wavelength of the fastest growing mode 
for non-axisymmetric MRI (left) and the wavelength divided 
by the grid resolution (right) on a meridional plane (x-z) at 
t — tmrg ~ 15.0 ms for the lowest-resolution run. 


« O.OIMq. The torus wind would significantly con¬ 
tribute to kilonova/macronova in BH-NS mergers. This 
point should be investigated systematically in a future 
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FIG. 5. Profiles of rest-mass density with velocity fields 
(top-left), plasma /3 (top-right), thermal component of spe¬ 
cific internal energy (bottom-left), and sum of the Maxwell 
and Reynolds stress (bottom-right) on the meridional plane 
at t — tmrg ~ 50.6 ms for the highest-resolution run. In all 
the panels, the black curves denote contours with ut = —0.98 
(dashed), —1.00 (solid), and —1.02 (dotted), respectively. In 
the bottom-right panel, the red (purple) colored region indi¬ 
cates that the energy is transferred outward (inward). 
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